Effective mass theory of monolayer 5-doping in the high-density limit 
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Monolayer 5-doped structures in silicon have attracted renewed interest with their recent incor- 
poration into atomic-scale device fabrication strategies as source and drain electrodes and in-plane 
gates. Modeling the physics of 5-doping at this scale proves challenging, however, due to the large 
computational overhead associated with ab initio and atomistic methods. Here, we develop an ana- 
lytical theory based on an effective mass approximation. We specifically consider the Si:P materials 
system, and the limit of high donor density, which has been the subject of recent experiments. In 
this case, metallic behavior including screening tends to smooth out the local disorder potential 
associated with random dopant placement. While smooth potentials may be difficult to incorporate 
into microscopic, single-electron analyses, the problem is easily treated in the effective mass theory 
by means of a jellium approximation for the ionic charge. We then go beyond the analytic model, 
incorporating exchange and correlation effects within a simple numerical model. We argue that such 
an approach is appropriate for describing realistic, high-density, highly disordered devices, provid- 
ing results comparable to density functional theory, but with greater intuitive appeal, and lower 
computational effort. We investigate valley coupling in these structures, finding that valley splitting 
in the low-lying F band grows much more quickly than the F-A band splitting at high densities. 
We also find that many-body exchange and correlation corrections affect the valley splitting more 
strongly than they affect the band splitting. 

PACS numbers: 73.22.-f, 85.35. Be, 73.21.Fg, 85.30.De 



I. INTRODUCTION 

As transistors continue to shrink in size, they ap- 
proach a limiting regime where the electrons are con- 
fined and controlled over atomic length scales. Silicon- 
based devices have shown particular promise in this re- 
gard. For example, electrons on individual donors or 
traps have been probed inside conventional metal on in- 
sulator field-effect transistors, 1-4 and have been proposed 
as qubits for quantum computing architectures. 5-11 In 
other experiments, quantum dots with fewer than ten 
deterministically-positioned donors have been tunnel- 
coupled to proximal leads, 12 and single electrons have 
been confined using electrostatic top-gates. 13 In sev- 
eral experiments, individual spins have also been 
measured. 1,2 ' 14 While such technologies will have appli- 
cations for conventional computing, much of the recent 
progress in this area has been spurred by the quest for 
spin-based qubits. 5 > 15 ~ 18 

In this work, we focus on devices formed of degen- 
erately doped phosphorus in silicon. In the labora- 
tory, the silicon is masked by hydrogen atoms, which 
are lithographically patterned using a scanning tunnel- 
ing microscope. 19-22 Phosphorus atoms from phosphinc 
gas are then incorporated into unmasked segments of the 
top (monatomic) layer of silicon. A self-limiting growth 



mechanism leads to rather uniform doping densities cor- 
responding to one substitutional donor for every four 
atomic sites in the 5-layer. 23 The resulting devices are 
fully epitaxial. 

Realistic theoretical models of Si:P disorder have 
proven challenging, with calculated band structures 
known to depend very sensitively on the disorder model. 
In density functional calculations, for example, the num- 
ber, the splitting, or even the existence of energy bands 
depends on the nature of the symmetries, the place- 
ment of the donors and the size of the computational 
unit cell. 24 In addition to disorder, fabrication geome- 
tries play an important role in device operation. Special- 
ized structures in typical devices may range in size from 
nanometers to microns, 12 causing technical challenges for 
any theoretical treatment. It is not currently feasible to 
treat large, disordered devices with atomistic accuracy; 
hybrid approaches, however, may provide viable, self- 
contained solutions for such large-scale problems. Multi- 
scale methods, in particular, hold great promise. Exam- 
ples include the merging of tight binding and local density 
techniques. 25 There are also computational advantages 
to eliminating disorder effects within the doping plane. 
In this case, translational symmetry can be restored by 
averaging. 24,26 

In this paper, we develop a coarse-grained theory of 
(5-doped Si:P devices, consistent with the effective mass 
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approximation. Effective mass theory (EMT) provides 
an efficient means for analyzing large, complex systems, 
like tunnel-coupled devices. Recently, such methods were 
applied to the problem of few-electron quantum dots. 12 
The main modification to the bulk EMT, required for 5- 
doping, involves the uniform shifting of energy bands, up 
or down. These shifts account for the quantum confine- 
ment in the (5-doping potential. A closely related system, 
whose band structure can be understood in terms of band 
shifting, is the inversion layer. 27 An important distinction 
between inversion layer and Si:P devices is that the lat- 
ter have very high carrier densities, which leads to the 
occupation of multiple bands. Density functional meth- 
ods confirm this picture of simple energy shifts in disor- 
dered Si:P geometries. 26,28 It is especially important to 
note that the parabolic shape of the bands and the corre- 
sponding curvature (the so-called inverse effective mass) 
are largely unaffected by the energy shifts. This suggests 
a straight-forward modification of the bulk effective mass 
theory to incorporate the confinement effects, which we 
describe in detail. 

We consider several approaches for obtaining a two- 
dimensional (2D)-EMT in Si:P. The simplest approach is 
empirical. In this case, the relevant parameters in the 
2D-EMT, which we can think of as inputs to the theory, 
are obtained from a more rigorous, ab initio band struc- 
ture calculation. When possible, the parameters are ob- 
tained directly from experiments. We go on to describe 
the physical features and phenomena that may be com- 
puted within the EMT. These include many-body effects, 
exchange and correlation, and valley splitting. The ques- 
tion of disorder and donor placement is not a concern for 
the EMT because of the coarse-grained nature of the the- 
ory. The characteristic length scale associated with the 
coarse-graining is the effective Bohr radius. For quarter- 
monolayer doping, this length scale encompasses many 
donors. A jellium approximation for describing the donor 
charge is therefore appropriate to our problem. Local 
variations of the jellium density could be introduced into 
this theoretical framework; we do not, however, consider 
such problems here. 

The theoretical approach we employ is similar to the 
numerical effective mass theory of Scolfaro, et al. 29 who 
consider a periodic array of thick low-density (5-layers. 
The larger separation between their donors required a 3D 
treatment and produced valley splittings of the A-band 
in the range of 20-46 meV, which is more consistent with 
individual donors than high-density (5-layers, where the 
A-band is essentially degenerate. 24 These increased with 
doping density and are inconsistent with ab initio results 
for high-density (5-layers. 24 Here, we treat the experimen- 
tally relevant case of thin layers of high doping density, 
where it is possible to project the 3D-EMT onto a 2D 
theory of immediate interest for 2D devices. As much 
as possible, we focus on analytical (rather than numeri- 
cal) methods, which allows us to identify the underlying 
physics of the (5-layers. For example, we obtain a den- 
sity scaling theory. Later, we obtain numerical solutions 



that provide theoretical parameters for the 2D theory. 
The present analysis is formulated in terms a single 5- 
layer; although our geometry is not periodic, multiple 
layers could be treated by a straightforward extension of 
our theory. Rodriguez- Vargas et al. 30 also solved a non- 
periodic double-layer system numerically, although their 
approach is semi-classical where ours is fully quantum 
mechanical. 

The paper is outlined as follows. In Sec. II we identify 
input parameters and develop a description for the shift- 
ing and filling of the bands within the EMT. In Sec. Ill, 
we clarify the main concepts of the shifted band model 
by deriving the (5-doping EMT from a bulk, 3D-EMT. 
This provides a setting to discuss the various compo- 
nents of the theory, which are not normally associated 
with EMT, but may be easily incorporated. These in- 
clude many-body interactions, valley splitting, exchange 
and correlation effects. The utility of the effective mass 
method is demonstrated in Sec. IIIB, where an analyti- 
cal, variational theory of (5-doping is derived. This leads 
naturally to a scaling theory for the quantum confine- 
ment lengths and the energies of the different conduction 
bands. We also perform a more rigorous, numerical anal- 
ysis of the (5-doping problem, to obtain an alternative set 
of EMT parameters, which we compare to results from 
more microscopic derivations in Sec. IV. We conclude in 
Sec. V. The two appendices provide further details on 

(A) the numerical methods employed in our work, and 

(B) our exchange-correlation analysis. 



II. TWO-DIMENSIONAL EFFECTIVE MASS 
THEORY 

In this section, we describe the modifications to the 
bulk conduction band structure of Si, due to (5-doping in 
the z — plane. For n-type devices in the low tempera- 
ture regime, the active electrons tend to fill only the low 
energy portions of the conduction band, known as val- 
leys. As is well known, 31 the valleys in bulk silicon are 
six-fold degenerate, with minima occurring in the equiv- 
alent [100] directions, about 85% of the way to the Bril- 
louin zone boundary. A given valley minimum therefore 
occurs at k ~ 0.85(27r/a), where a = 0.543 nm is the 
length of the cubic unit cell. To a good approximation, 
the low energy band structure in a given valley appears 
parabolic. For example, for the +x valley (along [100]), 



E +x ~ 



h 2 (k x - k ) 
2mi 



h 2 k 2 y h 2 k 2 z „ 



2m t 



2m 



(1) 



where m; ~ 0.92mo is the longitudinal effective mass, 
m t — 0.19too is the transverse effective mass, m is the 
bare electron mass, and E c is the conduction band mini- 
mum. A cut through E +X (k) along [100] is sketched as a 
dashed, red curve in Fig. 1(a). Here, we have set E c = 0, 
defining the band minimum as the zero of the energy. 
Although the band structure associated with the y and z 
valleys lies outside the range of the plot in Fig. 1(a), their 
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wavevectors also have k x components, which can be pro- 
jected onto the k x axis as shown with a dashed, blue line. 
These projected valleys are centered at k x — 0, analogous 
to the transverse terms in Eq. (1). The anisotropic ef- 
fective masses, mi and m t , together with the valley min- 
ima at fco, capture the main low-energy physics of this 
problem, and they form the main inputs to the 3D bulk 
effective mass theory. 

A single n-type donor ion, such as P, creates a local 
dip in the electrostatic potential, with a corresponding 
low-energy bound state 46 meV below the conduction 
band. For (5-doping in the z — plane, an electronic 
wave function extends over many donors. The electro- 
static potential and corresponding binding energies are 
much deeper than for a single donor. Because the elec- 
tron covers many randomly placed donors, it is conve- 
nient to ignore their individual positions, and to instead 
treat the dopants through a 2D "jellium" approximation, 
corresponding to a totally uniform charge distribution in 
the z = plane equal to the average 2D charge density of 
the discrete dopants. For an infinite sheet of charge, the 
electrostatic potential does not vary in the lateral plane. 
Figure 1(b) shows a vertical cut V(z) through such an 
electrostatic potential, together with the resulting con- 
fined states. Assuming overall charge neutrality, V(z) 
flattens out far away from the doping plane. The correct 
binding energy is obtained by aligning this asymptotic 
potential with the bottom of the bulk conduction band, 
as shown in the figure. 

The result of the vertical confinement V(z) is to re- 
duce the continuum of bulk Bloch states from 3D to 2D. 
Specifically, the allowable k z components of the Bloch 
states become quantized to form bound states. The 
parabolic energy decomposition of Eq. (1) suggests that 
the third term on the right-hand-side (the z term) should 
be replaced by a quantized band energy: 



E. 
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E A . 



2m; 2mt 
Similarly for the +y and +z valleys, we have 
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Analogous equations are obtained for the —x, —y, and 
— z valleys by replacing k —k . The quantized ener- 
gies E-p and E A depend on the effective masses in the k z 
terms of the bulk equations. The x-y valleys are 8-fold 
degenerate, including spin and valley degeneracies, while 
the z valleys are 4-fold degenerate. Later, we will discuss 
the lifting of the T degeneracy, in a process known as val- 
ley splitting. (For most cases of interest, the lifting of the 
A degeneracy, if present, will be negligible.) This leads to 
the distinct quantized energies, Ei and E 2 , correspond- 
ing to the Ti and T2 bands, which are shown in Fig. 1(b). 
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FIG. 1: (Color online) (a) Effective mass theory for bulk Si 
(dashed lines) and for 5-doped Si:P (solid lines). The minima 
of the bulk conduction bands define the energy zero, (b) The 
donors in the z — plane produce a laterally-averaged elec- 
trostatic confinement potential, with the predominant eigen- 
states Ti, F2, and A. The confinement along z leads to trans- 
verse bands along k x and k y that are shifted downward into 
the gap, as shown in (a). For all the plots, we take k y — 0. For 
the bulk r bands (dashed blue curve) we also take k z — ko, 
while for the bulk A bands (dashed red curve) we take k z — 0. 
The various bands are then projected onto the same axis. For 
the (5-doped bands, the energy shifts are determined by the 
eigenstates shown in (b). 



The leading order effect of vertical confinement is there- 
fore to reduce the 3D band structure to 2D, as shown in 
panel (a), with the 2D bands shifted downward by their 
respective binding energies. 

In the arguments presented above, the effective masses 
in the 2D theory should be identical to the bulk effective 
masses. Any deviations from the energy decomposition 
of Eq. (1) would imply mixing the effective masses and 
weaken the theory. Rigorous ab initio band structure 
calculations with a laterally averaged charge distribution 
confirm that the 2D effective masses for Si:P with the 
most prevalent 1/4 monolayer (ML) doping level are al- 
most identical to the bulk masses in Si, 26 and that the 
bulk bands are simply shifted downward by their respec- 
tive binding energies. 28 In short, all the evidence suggests 
that 2D and 3D effective mass theories should both be 
accurate in this system. 

The main parameters characterizing the 2D effective 
mass theory are Ei, E 2 , E&, m t , mi, and fco- Prefer- 
ably, their values should be obtained from experiments, 
or if not, from accurate ab initio theories. In a later 
section, we will show that E\, E 2 , and E& can be de- 
rived directly from the 3D effective mass theory. Addi- 
tional derived quantities of interest include the fractional 
fillings of the different conduction bands. These fillings 
have important implications for processes like transport 
and tunneling, which may occur much more readily in 
one band than another. While the total filling of the 
2DEG is determined by the total number of electrons, 
or the ionic charge density (assuming charge neutrality), 
the problem of calculating fractional band fillings is much 
subtler, since it depends on the accuracy of the binding 
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energy calculations. 

The respective filling fractions are defined as 0i, f3 2 
and (3 a, where 



(3l+(3 2 +(3 A = 1. 



(5) 



Alternatively, combining the Ti and T 2 bands into a sin- 
gle T band gives 



/3r+/?A = 1. 



(6) 



The corresponding charge densities are given by <r 7 = 
— /3 7 cr, where 7 = 1, 2 or A, and a is the average 
ionic charge density. For 1/4 ML doping, we have 
a = 0.2717 C/m 2 . Conventional 2D band filling argu- 
ments lead to 



Ep — E2 — 



Ep — E\ = (3a 
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(7) 



(8) 
(9) 



where we assume a shared chemical potential Ep at zero 
temperature. (Non-zero temperatures could also be con- 
sidered, although very low temperatures are assumed 
here, as appropriate for the main applications of inter- 
est.) For a combined T band we have 



Ep — Ey = 



irah 2 
2em t 



(10) 



The linear system of equations (5), and (7)-(9) [or (6), 
(9), and (10)] may readily be solved to obtain the filling 
fractions and Ep. For example, the T-A solution is given 

by 



E A + Er^mt/Ami + (irati 2 /Ae^m t mi) 
&f = ; ; , (11) 



/5r 



1 + T s /m t /Am l 
Ea — E r + (irati 2 /Ae^/m t mi) 
1 + y/mt/Ami^ (irati 2 /2em t ) 



(12) 



with p A = 1 - Pr- 

To conclude this section, we note that the 2D effective 
mass theory, described above, includes only the six low- 
lying valleys of the bulk conduction band structure. It 
is known that other bands may also begin to fill at the 
1/4 ML doping level; most notably, the IX /2X band may 
dip slightly below the Fermi energy. 26 These bands could 
be included in our 2D theory using the same methods 
described above. Their fractional fillings, however, are 
small enough that the main physics is already captured 
in the theory as presented here. 



III. THREE-DIMENSIONAL TREATMENT OF 
(5-DOPING 

In this section, we use the 3D-EMT to study the prob- 
lem of (5-doping in the z — plane. Due to our assump- 



tion of uniform doping in the lateral plane, the calcula- 
tion is one-dimensional in the variable z. Our ultimate 
goal is to derive the input parameters for a 2D-EMT. 

The theory must include many-body effects, and we be- 
gin by developing a simple Hartree theory. We go on to 
obtain a variational solution to the problem of (5-doping 
in Si, as well as scaling estimates for the vertical con- 
finement lengths and energies in the T and A bands. We 
show how valley splitting can readily be included as a cor- 
rection to the effective mass theory. Finally, we extend 
the theory to include exchange and correlation contribu- 
tions, which are used to obtain more accurate results for 
the 2D-EMT. 



A. Hartree Theory 

Within the jellium approximation, the ionic charge 
density is given by 



Pi(z) = ad(z), 



(13) 



where we take z = as the (5-doping plane. In an effective 
mass-Hartree theory, 31 the electron charge densities in 
the T and A bands are defined as 



Pr(z) = 
Pa(z) = 



-afcF 2 (z), 
-a0 A Fl(z), 



(14) 
(15) 



where Fr(z) and Fa(z) are the corresponding envelope 
functions obtained by solving the Schrodinger-like equa- 
tions 



ti^^_ 
2mi dz 2 
ti 2 d 2 
2m t dz 2 



+ V(z) 
+ V(z) 



F T (z)=e T F T (z), (16) 
F A (z) = e A F A (z). (17) 



Note that the full electronic wavefunctions also involve 
fast oscillations (e.g., Bloch oscillations) that occur over 
atomic length scales. These fast oscillations do not en- 
ter the envelope function equations (16) and (17). 32 In 
Sec. HID, below, we investigate perturbative corrections 
to the energy that occur when fast oscillations are taken 
into account in the T bands. At the present level of ap- 
proximation, however, Ti and T 2 have the same effective 
mass and the same envelopes. For now, we therefore con- 
sider just a single T band. Equations (6) and (13)-(15) ex- 
plicitly satisfy charge neutrality when the envelope func- 
tions are properly normalized. The full 3D charge density 
is then given by p(z) = pi(z) + p T (z) + p A (z); by design, 
we obtain J p(z) = 0. 

The electrostatic potential is calculated from Poisson's 
equation and Eqs. (13)-(15). It is convenient to compute 
the electric field contributions from each of the different 
charge sources. Making use of charge uniformity in the 
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lateral plane, we obtain 



Variational Calculation 



Ei(s) 
E 7 (z) 



Z(T 

Ye 
z 



sign(z), 



/ p 7 dz, 



(18) 
(19) 



where we have adopted the convention that / °° S(z) dz — 
1/2. Note that we will adopt the dielectric value e = 
11.4 eo throughout this work, as appropriate for silicon 
at very low temperatures. 

The corresponding contributions to the electrostatic 
confinement potential are then given by 



Vi(z) 



V 7 (z) 



-f 

Jo 



ea\z\ 
2e : 



E lz {z) dz, 



(20) 
(21) 



where E~. 



E • z. The total Hartree potential, used 



in Eqs. (16) and (17), is then given by V(z) = Vi(z) + 
Vr(z) + V A (z). Note that we have adopted the energy 
normalization V(0) — Vi(0) = 0. In this way, the po- 
tential minimum remains anchored, and is not affected 
by the specific details of the electronic wavefunctions. 
Such considerations simplify our variational calculation 
in Sec. IIIB, and are analogous to the Fang-Howard pro- 
cedure used for inversion layers. 31 (Later, for presenta- 
tion purposes, we will renormalize the energy as in Fig. 
1 such that the electrostatic potential is aligned with the 
bulk conduction band, in the region far from the 5-doped 
layer.) In principle, the methods described here could be 
modified to include external gates - for example, by in- 
troducing an external electric field. Wc do not consider 
this possibility here. 

In the Hartree many-body method, wc must obtain 
self-consistent solutions for the envelope functions (16) 
and (17), the Hartree potentials (20) and (21), and the 
Fermi level (11). The electron energies Ep and E A that 
appear in Eqs. (11) and (12), however, correspond to 
band minima; they are not the same as the single-particle 
energies £r and e A that appear in Eqs. (16) and (17). The 
band minima are given by 

E T = (T[m ; ]) r + (y j ) r + i(T/ r )r + (V A )r, (22) 
E A = (T[m t ]) A + (V,) A + (V r ) A + ^(V A ) A , (23) 

where T[m*] are the same kinetic energy operators ap- 
pearing in Eqs. (16) and (17), and the subscripts T and 
A refer to single-particle wavefunctions used to com- 
pute the expectation values. The prefactors of 1/2 in 
the Hartree terms prevent over-counting of the electron- 
electron interactions. 31 We emphasize that Eqs. (22) and 
(23) describe the band minima, and do not include any 
lateral kinetic energy. The lateral kinetic energy appears 
explicitly in the Fermi level equations. 



One of the main benefits of an effective mass theory 
is its simplicity and the ease with which solutions can 
be obtained. We take advantage of this now to obtain 
initial estimates for the band minima and the electron 
wavefunctions, by means of a variational method. We 
may even obtain simple analytical estimates, which allow 
us to scale the electron eigenf unctions and energy values. 

The problem can be formulated in several ways. Here, 
we consider a simple variational form for the single- 
electron wavefunctions, which is generally found to be 
consistent with more accurate treatments: 



F,(z) = 



( 2 



1/4 



(24) 



The wavefunction widths ar and a A , and the filling frac- 
tions pY and (3 Al represent the variational parameters in 
this approach, although we will use Eq. (6) to eliminate 
one of these variables (/3a)- The Gaussian form is partic- 
ularly effective in such a variation calculation because of 
its simplicity, and because it captures the essential prop- 
erties of the wavefunction (the width), while the maxi- 
mum value of the wavefunction is correctly determined 
via normalization. The Gaussian tail decays too quickly 
compared with a more realistic wavefunction; the tail, 
however, contributes very little to the expectation values 
in Eqs. (22) and (23), and therefore does not affect the 
leading order results of the variational calculation. 

Equation (24) immediately leads to analytical forms 
for quantities of interest, including the confinement po- 
tentials, 



V,(z) 



((Ta - ! -r (e-^K -l) 

I s| erf (\/2>|/a 7 ) , (25) 



eV87r 



e<7/3 7 . 



2c 



where erf(:r) is the error function. Equations (22) and 
(23) then reduce to 



Er 



2m; dp eV8n 
ea 



(l-Pr)yj4 + a 



(26) 
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o 7=ft\/ar + a 2 A 



(27) 
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3-72 y/2-1 \ . 
— ^ 1 ^ — Pr oa + Prar 



in terms of the variational parameters. 

To complete the variational analysis, we must mini- 
mize the average electron energy with respect to the vari- 
ational parameters. Here, we take the slightly different 
approach of minimizing the band energies Er and E A , 
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while introducing a level filling constraint from Eq. (12). 
The constrained problem is then converted to an uncon- 
strained problem by employing a Lagrange multiplier A. 
The minimization statement becomes 



Parameter value Eq. # 



V(E r + \g) = 0, 



(28) 



where the derivative is taken with respect to variables 
dr, oa, Pr, and A, and we have defined 



9 — E A -E r - 



■nah 2 



~Pr 



mt \ irah 2 

— - ■ (29 
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Note that the A derivative is equivalent to setting g = 0. 

Eliminating the Lagrange multiplier from Eq. (28) 
leads to a system of three equations 



dar 



dE A 
dPr 



dE r dE A 
dar da A 
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Ami 
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dE r dE A 
da A dar ' 
nah 2 



2em t 



dE r dE A 
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(30) 
(31) 

(32) 



which may be solved to obtain estimates for the varia- 
tional parameters. We do not report on such an analysis 
(yet), since Eqs. (30)-(32) cannot be solved exactly by 
analytical methods, and since we will perform a more 
rigorous numerical analysis later, which includes other 
contributions to the physics. The variational formula- 
tion, however, leads immediately to an important scaling 
theory, which we discuss now. 



C. Scaling Theory 

Our simple variational theory describes the main por- 
tion of the wavefunction envelopes correctly, and should 
therefore capture the leading order physics of the 5- 
doping problem. Based on this statement, we may draw 
some very general conclusions, which can be expressed in 
terms of a scaling theory. Of particular interest, the scal- 
ing theory captures the principal dependence of various 
quantities of interest, regarding the doping density a. 

The main expressions entering the variational proce- 
dure are given in Eqs. (26), (27), and (30) [or (12)]. We 
can reformulate Eqs. (26) and (27) in terms of dimen- 
sionlcss variables as follows: 



h 2 eV8-K 



1/3 



\ em a I 
( e 2 a 2 H 2 \ 1/3 



(33) 
(34) 



where the quantities with tildes are dimensionlcss, and 
mo can be taken as the bare electron mass. Recall here 
that E~ refers to the ground state confinement energy of 
the r or A band, as measured from the bottom of the 
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1.23 
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2.49 
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Er 


1.34 


34 


E A 


2.13 


34 
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0.13 


36 
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0.89 


37 



TABLE I: Dimensionless parameters appearing in the scaling 
theory. All scaling parameters were determined numerically, 
from the case of 1/4 ML doping, as described in Sec. IV. 



confinement potential. When the energy is normalized in 
this way, E 1 is strictly positive. 

Since (mt/mo), (to; /mo) and Pr are all of order unity, 
we expect that a 7 and E 7 should also be of order unity. 
Indeed, we may go beyond the variational calculation 
described above to obtain more rigorous numerical esti- 
mates (described below), for the case of 1/4 ML hlling. 
The results are shown in Table I. Note that these esti- 
mates could be improved by using results from rigorous 
microscopic calculations, or from experiments. The scal- 
ing theory itself, however, would remain unaffected. 

Based on Eqs. (33) and (34), we can deduce the scaling 
behaviors for other quantities of interest. For example, 
from Eq. (12) we obtain the relative filling fractions 



Pr ^ 0.19 + 6 



ah 4 ir 4 e 2 



Pa - 0.81 - I ( 



\ah 4 ir 4 



1/3 



1/3 



(35) 
(36) 



In the Hartree theory, the electrostatic potential is de- 
fined as V(z) = Vi(z) + V r (z) + V A (z). The depth of 
the confinement potential, Vo = V(oo) — V(0), plays an 
important role for the quantum theory. Within the vari- 
ational approach described above, this quantity can be 
expressed as 



V 



ecr 



ev/87r 



{a r Pr + a A p A ) . 



(37) 



To conclude this section, we note that the scaling the- 
ory breaks down outside a regime of validity. In the 
present analysis, we have assumed a jellium model for 
the doping. In the low density limit, however, the jel- 
lium m odel b reaks down when the average dopant separa- 
tion yje/ira approaches the characteristic effective mass 
length scale, min[ar,aA]- Within the scaling theory, we 
can estimate this breakdown density as 1/18 ML. In the 
high density limit, it is important to note that we have 
only included the F and A bands in the present analy- 
sis. For densities larger than 1/4 ML, the filling of ad- 
ditional bands becomes important. This can easily be 
accomplished and incorporated into the present formal- 
ism although it lies outside the scope of the present work. 
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D. Valley Splitting 

As discussed in Sec. II, the combination of inho- 
mogcneous (vertical) confinement and effective mass 
anisotropy lifts the degeneracy of the bulk bands. The 
resulting splittings can be quite large. Remaining degen- 
eracies are lifted when the confinement potential is sharp. 
For patterned devices, the valley splitting in the A band 
is extremely small, due to weak lateral confinement. 12 
In this section, we focus on the coupling between the z 
valleys, due to the sharp 5-doping potential, which splits 
the r band to form Ti and T2 bands. 

The envelope function equations (16) and (17) do not 
explicitly take into account the fact that the envelopes 
are formed from Bloch states within a given valley. We 
can account for this translation in the Brillouin zone in a 
simple way, by introducing an overall phase factor. 32 In 
the absence of valley coupling, we may therefore define 
the z- valley basis as follows: 33-35 



f ± (z) = e ±ik °*F r (z). 



(38) 



Here, Fy(z) is obtained from Eq. (16), and the resulting 
basis states f±(z) are degenerate. Valley coupling can 
be treated perturbatively in the same basis, through the 
Hamiltonian 



vo 



Vvo 



(39) 



There are two types of contributions which enter the 
valley-orbit coupling term Vyo- The first type can- 
not be described within an EMT. These include the so- 
called central cell corrections, which arise due to core 
electrons, 36 as well as discretization effects associated 
with the crystal lattice. 37 The latter contributions are 
fairly weak. Central cell corrections are also weak for 
shallow donors such as Si:P. Indeed, for isolated donors, 
the main contributions to valley-orbit coupling may be 
treated effectively using methods similar to those de- 
scribed below. 38 For simplicity, we therefore ignore non- 
EMT corrections here. 

Instead, we focus on valley-orbit couplings, which may 
be treated perturbatively within the EMT. They are de- 
fined as 



Vyo = (+\V\-) 



r 

j — ( 



V{z)F^{z)e- 2lkaZ dz 1 (40) 



where V(z) is the electrostatic confinement potential. 
The valley-split single-electron energy levels are then 
given by £1 — £r — I Vyo I and £2 = £r + |Vvo|i while 
the valley splitting is given by 2|Wo|- Similarly, the in- 
dividual band minima are given by 



E 1 =E T -\V vo \, E 2 = E T + \Vvo\- 



(41) 



The self-consistent numerical solutions, described be- 
low, arc unaffected by valley splitting. This becomes 



clear if we note that while the Ti and T2 bands fill dif- 
ferently, the electrostatic equations depend only on the 
combined filling factor, 0r = /?i+/?2- Likewise, the quan- 
tum mechanical envelope function equations are identical 
for Ti and T2- We may therefore solve for the energies 
Ey and and the fillings /?r and /?a, as we did pre- 
viously, while computing the perturbations due to valley 
splitting a posteriori. After solving for the envelope func- 
tions, Vyo is determined from Eq. (40). The band min- 
ima are then obtained from Eq. (41), while the fillings 
are obtained from 



/3 r , em t \Vvo\ 

31,2 = -7T ± 



(42) 



It is interesting to analyze the scaling behavior of 
the valley splitting, since energy splittings can be mea- 
sured experimentally, via spectroscopy techniques. Since 
ko ~ 1/a, Eq. (40) may be regarded as an integral trans- 
form that picks out the Fourier components in V(z)F 2 (z) 
with very short wavelengths. The predominant feature 
at short wavelengths is the sharply peaked confinement 
potential at z = 0. The more slowly varying features 
occurring away from z — effectively cancel out, and do 
not contribute greatly to the integral. (In some cases, a 
sharp variation of the wavefunction can also contribute 
to Vyo j 39 however, we ignore this possibility in the sim- 
ple estimate presented here.) We may therefore ap- 
proximate Vyo by truncating the integration range in 
Eq. (40) to a single oscillation of the exponential, from 
z = — 7r/2fc to 7r/2fc . Over this range, we can approxi- 
mate V{z)Fy{z) ~ ecr\z\/\/2irear, leading to the follow- 
ing estimate for the valley splitting: 



2|V vo | 



2tt ear 



(43) 



Despite the obvious roughness of this estimate, we expect 
it to encompass the leading order contributions to the 
scaling theory, which we find to be 



2|V vo | 



- ( m ° 



4 4 
e a 



V87r 2 fi 2 e 4 fc 6 



1/3 



V. 



(44) 



Thus, we note that the ri-^ splitting exhibits a much 
stronger dependence on the doping density a than does 
the IVA splitting, as borne out by the numerical analy- 
sis, described below. 



E. Exchange and Correlation 

Self-consistent many-body effects were included in the 
variational calculations of Sec. IIIB. There, we employed 
a Hartrcc theory, in order to simplify our analytical cal- 
culations. For more accurate numerical results, it is im- 
portant to also include exchange (X) and correlation (C) 
effects. We do this here, using the local density approxi- 
mation (LDA) . 40 Although we are primarily interested in 
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FIG. 2: Electronic charge density for the case of 1/4 ML 
doping in a Si:P r5-layer. Black solid line: density functional 
theory from Ref. 28. The density functional results exhibit 
small oscillations, which arise due to perfect ordering in the 
2D dopant array. Gray line: variational calculation, from 
Sec. IIIB. Dashed black line: self-consistent numerical so- 
lutions for the wavefunction envelope, as described in Ap- 
pendix A. Note that we have assumed a smooth lateral dis- 
tribution for the dopants, via the jellium model. The resulting 
effective mass solutions do not oscillate. Also note that the 
Gaussian variational solution accurately represents the full 
numerical results, except in the tail region. 



higher densities, it is well known that exchange and cor- 
relation effects are most significant at lower densities. For 
completeness, we will therefore study a wide range of den- 
sities, over which the LDA must remain valid. Typically, 
this requires specially developed parameterizations. 41-43 
Here, we follow the method of Scolfaro et al. 29 and 
Rodriguez- Vargas & Gaggero-Sager, 30 where the LDA is 
applied to the EMT equations, rather than at the atom- 
istic level, although we use a more recent and accurate 
density functional parameterization developed by Perdew 
and Wang 43 . The final outcome is a pair of new poten- 
tial terms, Vx(z) and Vc(z), which we add to the total 
electron confinement potential V(z): 



V =V + V X + V C . 



(45) 



lation potentials may be approximated as 

d(ne x ) _ 



Vx(z) = 
V c (z) = 



dn 
d(ne c ) 
dn 

-2A 



(4neh) 2 



9n\ 1/3 1 



nr. 



{Anehf 



{( 



1 + ^lln 



1 + 



(46) 
(47) 

1 



2Af 



(l + air a ) /' 



2Af)J' 



where the 3D particle density is given by n(z) — p(z)/e, 
and we define 



r s (n) 



Ana* 3 n 



-1/3 



f(r s ) = b^J 2 + b 2 r s + b 3 r 3 / 2 + 6 4 r s 2 . 



(48) 
(49) 



Following Perdew and Wang, the exchange and corre- 



Following Refs. 29 and 30, we adopt a geometrically av- 
eraged effective mass, m* = (m^m;) 1 / 3 , and an averaged 
effective Bohr radius, a* — (Aneh 2 /m*e 2 ). The param- 
eterization constants we use are appropriate in the ab- 
sence of spin polarization: A = 0.031091, ct\ = 0.21370, 
&i = 7.5957, b 2 = 3.5876, b 3 = 1.6382, and b 4 = 0.49294. 

We can estimate the magnitude of the exchange and 
correlation terms. We specifically consider the case of 
1/4 ML <5-doping. Because of the high doping density, 
we find that the dimensionless electron separation length 
r s ranges from about 0.6 at the center of the <5-doping 
layer to oo far away from the doping plane. In the high 
density region, which is our main interest here, the ex- 
change potential dominates over the correlation poten- 
tial. We can estimate its depth directly from Eq. (46), 
finding that Vxo — 70 meV. This may be compared to 
the electrostatic potential depth in Eq. (37), which we 
find from numerical calculations to be Vb — 670 meV. 

In the numerical calculations discussed below, we solve 
the 5-doping problem using the 3D EMT, including val- 
ley splitting, exchange and correlation effects, as outlined 
in Appendix A. Overall, we find that exchange and cor- 
relation have relatively small effects on the population of 
the bands or on the valley-splitting. The z-confinement 
of the electrons, however, is enhanced, particularly in 
the A-band. This is to be expected, since the exchange 
and correlation interactions both deepen the potential 
well, particularly near the doping layer, where the elec- 
tron density is highest. Since the A wave functions still 
spread out well beyond the doping plane, the overall ef- 
fect of exchange and correlation on the T-A splitting is 
fairly weak. The effect on the Ti-r2 valley splitting is 
stronger, however, as it is directly related to the sharp- 
ness of the confinement potential. Additional discussion 
of the exchange and correlation contributions to <5-doping 
is presented in Appendix B. 
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IV. RESULTS AND COMPARISON 

In this section, we describe the numerical results of 
our 3D-EMT description of ^-doping in Si:P, with details 
given in the Appendices. We also compare our results to 
other reports in the literature, mainly based on density 
functional theory. The closest points of comparisons in- 
clude the planar Wannier orbital (PWO) method of Qian 
et al., 26 and the full density functional calculations of 
Carter et al., 24 using the "single-zeta plus polarization" 
(SZP) basis set. Both of these techniques may be mod- 
ified to include disorder effects. For the PWO method, 
this was accomplished using laterally averaged confine- 
ment potentials. For the SZP method, quasi-random 
dopant arrays were considered, as well as "mixed" pseu- 
dopotentials, averaged over the atoms in the doping layer. 
In the interests of brevity and generality, we compare ex- 
plicitly to Ref. 26, and also to the full ab initio results 
of Ref. 28. Other abbreviations used hence are as fol- 
lows: our effective mass theory (EMT), the PWO method 
including short-ranged interactions between the dopants 
(PWOf), 26 and the fully ordered (SZPo) vs. partially 
disordered dopant arrays (SZPd) discussed in Ref. 28. 

Before discussing our main numerical results, it is im- 
portant to emphasize that the parabolic band structure 
assumed in the EMT (e.g., Fig. 1) is far more consis- 
tent with the highly disordered implementations of the 
density functional theory. The small unit cells associ- 
ated with dopant ordering generate effective terms in the 
Hamiltonian which couple the donor bands and lead to 
band structures that differ greatly from bulk silicon. 28 
We therefore conjecture that the EMT based on the jel- 
lium donor model should be understood as a highly dis- 
ordered model. Dopant ordering, or any specific type of 
disorder, can be introduced into the EMT through ad- 
ditional modifications of the jellium model. The smooth 
charge profiles obtained by EMT in Fig. 2 are also consis- 
tent with spatial averaging in the presence of disorder, as 
compared with the more oscillatory profile obtained from 
density functional theory for an ordered dopant array. 

For the case of 1/4-ML doping, we can summarize 
our main EMT results as follows. A fit of the numer- 
ical wavefunction envelopes to the Gaussian form used 
in the variational procedure of Sec. Ill B gives the enve- 
lope widths ar = 0.64 nm and a a — 1-30 nm. The band 
filling parameters are given by /?i = 0.19, /?2 = 0.18, 
and /3a = 0.63, while the valley splitting is given by 
2|Vyo| = 19 meV. In Figs. 3-5, we plot our numerical re- 
sults for the band minima and the wavefunction 
widths (in terms of the Gaussian fitting parameters ar 
and a a), and the band filling fractions /3r and /3a, as a 
function of doping density. The figure insets provide ad- 
ditional comparisons with the literature. (In some cases, 
the values have been determined graphically, from pub- 
lished plots.) 

We now discuss the main plots in Figs. 3-5 in more 
detail. In each plot, the markers represent numerical 
results obtained as a function of P doping density in the 
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FIG. 3: (Color online) Minimum band energies relative to the 
bulk band minima for the Ti band (blue circles), P2 band (red 
squares) and A band (green diamonds). Scaling theory fits are 
also shown as dashed lines. The minima of the confinement 
potential, Vo, are shown as black crosses. The corresponding 
scaling theory for the confinement minima Vo', based on de- 
rived parameters, is shown as a dashed black line (see main 
text). A rescaled fitting is shown as a dashed grey line. Inset 
shows a comparison of 1/4 ML results with Refs. 28 (SZPd) 
and 26 (PWO, PWOf) for the minima of the Ti (triangles), 
T2 (asterisks), and A (crosses) bands. (See main text for 
abbreviations.) 



6 layer, while the curves reflect the corresponding scaling 
theory parameters given in Table I. Direct comparisons 
are made to Ref. 26 (a fully-disordered technique), and 
to Ref. 28 (a full ab initio technique). 

Figure 3 shows the energies of the various band min- 
ima measured from the bottom of the bulk (undoped) 
Si conduction band. In these calculations, we have not 
considered background dopants, so the bulk band mini- 
mum (V — 0) corresponds to the asymptotic value of the 
confinement potential V(z) far from the doping layer. As 
the doping increases, V(z) deepens significantly, dragging 
the confinement energy levels with it. The valley split- 
ting between the T\ and T2 band minima also increases 
quickly at higher densities. The inset shows agreement 
between EMT and the SZPd method for the A and T 2 
bands, and agreement with the PWO and PWOf meth- 
ods for the valley splitting, for the 1/4 monolayer doping 
density case. 

We observe that the scaling theory describes the loca- 
tion of the band minima quite well. (The scaling of the 
valley splitting will be discussed below.) It appears that 
the scaling theory is less accurate for the potential well 
minimum Vq = V^(oo) — 1^(0). In this case, the scaling 
form (dashed black line) was derived from Eqs. (33)-(37) 
using the parameters in Table I. This discrepancy in the 
scaling theory is primarily due to the inclusion of cor- 
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FIG. 4: (Color online) Gaussian widths of the wavefunction 
envelopes: ar (blue circles), aA (red squares). Scaling the- 
ory fits are shown as dashed curves. Inset shows comparison 
of results for the full- width at half-maximum of the full elec- 
tronic density: this work (crosses), SZPo results from Ref. 28 
(asterisks). 



relation and exchange effects in the numerical solution, 
while they are absent from the discussion leading to the 
scaling theory. Exchange and correlation both tend to 
deepen the confinement potential, and they both have 
their greatest (absolute) effect at the origin. It is a sign 
of robustness of the scaling theory that such corrections 
can be accommodated by a simple adjustment of the scal- 
ing parameters, as indicated by the dotted black line. 

The a 7 parameters, or Gaussian widths for the wave- 
functions, are shown in Fig. 4. As expected, higher dop- 
ing tends to enhance the electron confinement, despite 
the Coulomb repulsion between the larger number of 
electrons. At low densities, the electron wavefunctions 
appear unphysically large, although our predictions for 
densities lower than 1/18 ML should not be compared 
directly to physical systems, as explained in Sec. IIIC, 
due to the breakdown of the jellium approximation. 

The scaling theory provides an excellent description of 
the numerical results over the entire range of densities 
in Fig. 4. Small inaccuracies of the scaling theory may 
be attributed to exchange and correlation effects, as dis- 
cussed in Appendix B. 

We expect the filling fractions plotted in Fig. 5 to 
match the predictions of scaling theory quite well, be- 
cause Eq. (12) is exact for a parabolic band structure, 
and because the band minima are also well described by 
the scaling theory. Note that we have plotted the results 
for fa and fa separately, and compared them to the scal- 
ing theory result for fa/2, where fa — fa+fa. As before, 
the main deviations from the scaling theory occur at low 
densities where exchange and correlation effects are most 
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FIG. 5: (Color online) Filling fractions for the different bands: 
fa (blue circles), fa (red squares), and /3 a (green diamonds). 
Scaling theory fits are shown as dashed curves (see main text) ; 
as discussed in Sec. HID, and in light of Eq. (42), half the 
scaling theory value for the F filling fraction (/3p/2) is plotted. 
Black, downward-pointing triangles correspond to /3a results 
from Ref. 26. Inset shows a comparison of fa (asterisks), fa 
(triangles), and /3a (crosses) for the case of 1/4 ML doping. 
(See main text for abbreviations.) 



important. 

Because Eq. (12) is generic, the scaling theory for fa 
could also be applied to results from other methods, such 
as those in Ref. 26. The latter (fa\ values) are graphically 
estimated and shown in the main panel of Fig. 5. The 
asymptotic, high density values of fa depend only on the 
effective mass, and should be nearly identical to those 
calculated here. The main deviations between the scaling 
theory and EMT at high densities arise because of small 
errors in the scaling theory for the quantity — Er). 
At low densities, the discrepancies are due to exchange 
and correlation effects. 

For the 1/4 ML results shown in the inset of Fig. 5 
(graphically estimated from bandstructures in the rele- 
vant papers) , our EMT results are most similar to PWOf. 
For the SZP results, the disordered model (SZPd) is most 
similar to EMT. This is consistent with our conjecture 
that the EMT provides a good description of the high 
disorder limit. We also note that, using the definitions of 
in Eqs. (7)-(9), it appears that the charge neutrality 
condition, Eq. (5), is not satisfied in either case (partic- 
ularly in Ref. 28) - though this is due, at least in part, 
to their inclusion of other bands, as discussed below. 

Figure 6 shows our calculated results for the Fermi 
energy, relative to the bulk conduction band minimum. 
Some Fermi energies obtained by density functional 
methods are also shown. Fermi levels are notoriously dif- 
ficult to calculate accurately. This is especially true for 
highly doped Si:P, due to the filling of multiple bands, 
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FIG. 6: Comparison of Fermi levels obtained in this work 
(circles), PWO results from Ref. 26 (diamonds), and SZPo 
results from Ref. 28 (squares). 
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FIG. 7: Numerical solutions for valley splitting between the 
Ti and bands (circles), with scaling theory fit (dashed line). 



and the fact that separate band minima must all be 
computed self-consistently. The EMT results in Fig. 6 
change sign, unphysically, near the 1/4 ML doping level. 
This can be attributed to the conspicuous absence of the 
IX /2X bands at this density, which we have chosen not 
to include in our model in light of the fact that the rele- 
vant filling fraction has been found to be less than 0.01 26 . 
Such high-lying bands would absorb high energy elec- 
trons, and would therefore lower the Fermi level as they 
begin to fill. Although no experimental measurements 
of the Fermi level are available at the present time, pre- 
liminary tunneling experiments between nano-fabricated 
wires suggest a Fermi energy of about -20 meV for the 
case of 1/4 ML doping. 44 We also direct the reader to 
further discussion in Appendix B. 

Finally, we plot our results for the valley splitting be- 
tween the Ti and T 2 band minima in Fig. 7. The valley 
splitting varies by several orders of magnitude over this 
density range. Carter et al. have noted that the valley 
splitting is particularly sensitive to the disorder model 
used in the calculations, with ordered dopants typically 
leading to larger valley splittings. 24 

The scaling results shown in Fig. 7 provide an excel- 
lent representation of our numerical solutions. It is in- 
teresting to note that, once again, the deviations are due 
to exchange and correlation effects (and partially to the 
truncation of the integral in Eq. 40 and subsequent lin- 
ear approximation to the potential in this region). In 
this case, however, the deviations are most evident at 
high densities. This occurs because exchange and corre- 
lation deepen the confinement potential at high densities, 
while the sharpness of the confinement potential provides 
the main contribution to the valley splitting. 



V. SUMMARY AND CONCLUSIONS 

In this paper, we developed an effective mass theory 
for high density (5-doped Si:P, and we argued that the 
model is consistent with the limit of high disorder. The 
method was applied to study infinite planes of Si:P. First 
a variational model was solved, which provided simple 
analytical results and demonstrated a remarkable agree- 
ment with density functional theories for very few as- 
sumptions. Second, a more comprehensive numerical 
model was solved, including exchange and correlation ef- 
fects, and valley splitting between the Ti and r 2 bands. 
Self-consistent solutions were obtained for systems com- 
prised of r5-layers with P densities ranging from 1/512 to 
1 monolayer. 

In our model the inclusion of valley-splitting in the 
self-consistent description has no effect on the extent of 
the donor wavefunctions, or on the relative populations 
of the r and A bands. Splitting between the Ti and r 2 
band minima naturally produces a difference in the filling 
of those bands. Similar to Ref. 26, we find that electrons 
are mainly concentrated in the T-bands at low densities, 
while nearly two thirds of the electrons shift to the A 
band for densities approaching 1/4 ML. We predict that 
the A filling should approach 70-80% at higher densities, 
although such densities are not easily achieved in physical 
systems. 

Since effective mass calculations mainly involve solving 
for coarse-grained envelope functions, they can be im- 
plemented much more efficiently than ab initio methods, 
such as density functional theory, or tight-binding. Speed 
comes at a price, however, since the technique relies upon 
accurate input data, including the anisotropic effective 
masses for the conduction electrons, the dielectric con- 
stant of the host material, and knowledge of the under- 
lying bulk band structure. The effective mass method 
may also have limited applicability at low and high den- 
sities. In the first case, the jellium approximation breaks 
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down below the metal-insulator transition. In the second 
case, the T and A bands tend to over-fill above 1/4 ML 
doping density, due to the absence of additional bands 
in the present theory. Fortunately, the applicable range 
includes most problems of current experimental interest. 
One of the main attractions of the effective mass theory 
is its versatility and its potential for treating complex 
problems of current interest for devices. 12 The method is 
easily extended to higher dimension and large-scale ge- 
ometries, which are typically out of reach for ab initio 
techniques. 

It is perhaps surprising that a minimal model like 
the effective mass theory could provide such a reason- 
able account of the broad range of physics in these 5- 
layers. The power of such a simple representation has 
been well-illustrated in the preceding sections, suggest- 
ing that the inputs to the effective mass theory capture 
the main physics in this problem. Our results are well 
matched to the predictions not only of Qian et a/., 26 
whose technique incorporates disorder in a similar man- 
ner as ours, but also to Carter et at, 28 who utilize quasi- 
disordered dopant arrays. Our calculated valley-splitting 
agrees with Rcf. 26, and our IVA band splitting and 
binding energies agree well with Ref. 28. Our results for 
the band fillings compare well to other values in the liter- 
ature; in particular to the fully-disordered results of Ref. 
26. 

Finally, we have condensed our effective mass results 
into a scaling theory, which may represent the simplest 
and most far-reaching outcome of the theory of infinite 
Si:P planes. The scaling theory reproduces our numeri- 
cal results very well up to high doping densities of order 
1/2 ML, and it enables analytical calculations of vari- 
ous physical quantities, as a function of the doping den- 
sity. Examples include the band energies, E 1 <~ er 2 / 3 , 
the Gaussian widths of the wavefunctions, a 7 <~ cr -1 / 3 , 
the depth of the confinement potential, Vq ~ a, and the 
valley splitting, 2Vyo ~ ct 4 / 3 . 
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Appendix A: ID Numerical Procedure 

In this appendix, we obtain full numerical eigenvalue 
solutions of Eqs. (16) and (17), in the presence of cor- 
relation and exchange. In contrast with the variational 
calculations of Sec. IIIB, we do not fix the value of (Hp- 
The method we use is iterative. At each stage of the 
procedure, a new density profile is obtained. The density 
is used to compute the electrostatic potential, which is 
used to compute a new density profile, and so on. Self- 
consistency is accomplished by introducing a new con- 
straint: 

h = J[n a+1 (z)-n a (z)] 2 dz = 0. (Al) 

Here, a and a + 1 indicate the iteration number. In other 
words, the density should not change from iteration to 
iteration once convergence is achieved. 

The self-consistency constraint of Eq. (Al) is well- 
defined and can be applied to density approximations 
involving many parameters. For example, the density 
could be defined spatially as = n(zi), where i is now a 
spatial index (not the iteration index, a). In this case, the 
self-consistency procedure is high-dimensional, involving 
many independent parameters. To simplify our anal- 
ysis, we will calculate the electrostatic potential using 
Gaussian density profiles, similar to those in Sec. IIIB. 
From Fig. 2, we see that such Gaussian forms provide a 
very reasonable estimate for the density, and can be im- 
mediately integrated to obtain Hartree potentials, as in 
Eq. (25). Indeed, by fitting Gaussian forms to "exact" re- 
sults for Fr and Fa, obtained by finite clement methods, 
we obtain density approximations that are much more 
accurate than the variational approximation shown in 
Fig. 2. 

Our self-consistent method is then expressed as follows: 
(i) provide a Gaussian estimate for the density profile in 
iteration a, (ii) incorporate the corresponding Hartree, 
exchange and correlation potentials into a finite element 
Schrodingcr solver, (iii) fit the resulting eigenfunctions 
F r and Fa to Gaussian forms, to be used in iteration 
a + 1, (iv) repeat until convergence is achieved, by mini- 
mizing the constraint of Eq. (Al), where n a refers to the 
Gaussian forms. The minimization step (i.e., the con- 
straint) is multidimensional in the Gaussian parameters 
ar and a a, and can be accomplished using the BFGS 
method. 45-48 Note that for a more accurate result, the 
wavefunction could be expanded in a larger basis of Gaus- 
sian functions. Minimization of such a parameterization 
would still be accomplished more efficiently than for the 
spatial parameterization, {n>i}. 

It is possible to incorporate the self-consistency con- 
straint of Eq. (Al) into the variational construction of 
Eq. (28). The quantity to be minimized would then be 
/ = Ep + Xg + fjbh. In such an approach, h = would 
never be satisfied until convergence is achieved. We have 
found, however, that the extended parameter space asso- 
ciated with allowing h ^ introduces new local minima, 
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which arc difficult to avoid. We therefore employ a dif- 
ferent method, as described below. 

The quantum mechanical problem involves two 
Schrodinger equations, (16) and (17), whose coupling, 
through the electrostatic potential V(z), is fully speci- 
fied by the parameter /?r- Consequently, the parameters 
a r and a a are completely determined by /3 r , and do not 
depend on parameters used in the variational approach, 
such as A. This statement remains true for more com- 
plex geometries, such as the 2D geometries considered 
in Ref. 12. In such cases, many parameters may be re- 
quired to fully describe the wavefunctions, although self- 
consistency still depends only on /3p. We can therefore 
use a numerical approach where the Schrodinger equa- 
tions are solved self-consistently for a fixed value of /?r- 
This one parameter is then varied, in order to satisfy the 
Fermi level constraint, g — 0. The latter problem is sim- 
ple, and can be accomplished using a Newton-Raphson 
method. 

The BFGS method, used to achieve self-consistency, 
involves solving for v different parameters, such as ar and 
a a , which may not have the same dimensions. The tech- 
nique involves calculating a v x v Hessian matrix whose 
elements may have values that differ by many orders of 
magnitude. It is numerically challenging to invert such 
a matrix. Therefore, to make the problem tractable, we 
transform to dimcnsionless variables. We have already 
identified the appropriate quantities for rescaling lengths 
in Eq. (33) and energies in Eq. (34). The envelope func- 
tions and Gaussian form for the electron density may be 
simply expressed in these terms, as may the potentials 
and the Fermi constraint. 

The Schrodinger equations are solved by finite element 
methods. The average dimensionless energy expecta- 
tion values Ey and Ea are used in the Fermi constraint. 
These are readily computed by our finite element solver, 
via Eqs. (16) and (17) as 

E r = cy-^(Vy)y, (A2) 

E a = e A - 1(Va)a- (A3) 

The parameters to be solved for in our ID model are 
then {ar,aA,/?r}- The valley splitting does not affect 
these solutions, and is therefore calculated post hoc, as 
described in Sec. HID. 

The final results are obtained numerically for the 1/4 
ML case, giving ar = 0.64 nm and oa = 1-30 nm. The 
value of /3r is slightly larger than the estimate /? r ~ 1/3 
obtained in Ref. 26, and larger than the estimated value 
of 0.310 from Ref. 28 (see Sec. III). 

It should be noted that, with appropriate initial 
guesses as to the input parameters, the model converges 
in a matter of minutes on a standard laptop. This is far 
more efficient than a typical ab initio calculation, which 
often requires tens of hours' runtime across tens of cpus, 
and could potentially be used to model much larger sys- 
tems, with device scales beyond the reach of more rigor- 



Dopant density 


ar 


a r 




OA 


(ML) 


(nm) 


(nm, no XC) (nm) 


(nm, no > 


1 


0.403 


0.412 


0.825 


0.876 


2/3 


0.461 


0.472 


0.943 


1.007 


1/2 


0.508 


0.521 


1.036 


1.112 


1/3 


0.581 


0.598 


1.183 


1.280 


1/4 


0.639 


0.660 


1.298 


1.415 


1/8 


0.804 


0.837 


1.625 


1.804 


1/16 


1.009 


1.063 


2.027 


2.302 


1/32 


1.265 


1.350 


2.521 


2.948 


1/64 


1.580 


1.718 


3.120 


3.774 


1/128 


1.959 


2.185 


3.854 


4.852 


1/256 


2.420 


2.788 


4.677 


6.208 


1/512 


2.955 


3.547 


5.601 


7.993 



TABLE II: Envelope widths (ar,A values) with and without 
exchange and correlation. 



ous techniques. 



Appendix B: Exchange and Correlation 

The purpose of this Appendix is to detail when inclu- 
sion of correlation and exchange effects in the numeri- 
cal model is necessary for accurate calculation. It as- 
sumes the treatment given above, accounting for valley- 
splitting. To characterize the effects of exchange and cor- 
relation, we perform our self-consistent calculations both 
with and without Vx and Vc in the Hamiltonian. 

The primary effect of including exchange and correla- 
tion (XC), as discussed in the main text, is to slightly 
deepen the potential well and steepen it in the vicinity 
of the minimum. This results in a contraction of the en- 
velope functions, as can be seen in Table II. The effect 
is more marked for the low-density cases, where the jel- 
lium approximation is in question and (Vxo + Vco) /Vq is 



Dopant 


Pi 


Pi 


P2 


P2 


Pa 


Pa 


density 




(no XC) 




(no XC) 




(no XC) 


(ML) 














1 


0.169 


0.169 


0.146 


0.147 


0.685 


0.684 


2/3 


0.174 


0.174 


0.157 


0.158 


0.670 


0.669 


1/2 


0.178 


0.178 


0.164 


0.165 


0.658 


0.657 


1/3 


0.186 


0.186 


0.175 


0.175 


0.640 


0.639 


1/4 


0.192 


0.192 


0.182 


0.183 


0.626 


0.625 


1/8 


0.210 


0.209 


0.202 


0.202 


0.588 


0.590 


1/16 


0.231 


0.228 


0.225 


0.223 


0.544 


0.549 


1/32 


0.255 


0.249 


0.251 


0.245 


0.494 


0.506 


1/64 


0.284 


0.272 


0.280 


0.269 


0.437 


0.459 


1/128 


0.317 


0.295 


0.314 


0.293 


0.370 


0.412 


1/256 


0.357 


0.319 


0.354 


0.317 


0.289 


0.365 


1/512 


0.408 


0.341 


0.406 


0.339 


0.186 


0.320 



TABLE III: Filling fraction (/3 values) with and without ex- 
change and correlation. 
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Dopant 


rvr 2 


rvr 2 


Ti-A 


ri-A 


density (meV) 


(meV) 


(meV) 


(meV) 


(ML) 




(no XC) 




(no XC) 


1 


183.2 


175.6 


595 


600 


2/3 


88.3 


84.0 


446 


449 


1/2 


54.1 


51.6 


363 


365 


1/3 


28.8 


27.6 


271 


272 


1/4 


19.0 


18.2 


220 


220 


1/8 


7.3 


7.0 


132 


132 


1/16 


2.8 


2.7 


79.2 


77.8 


1/32 


1.1 


1.0 


47.1 


45.4 


1/64 


0.44 


0.41 


27.9 


26.2 


1/128 


0.18 


0.16 


16.4 


14.8 


1/256 


0.07 


0.06 


9.7 


8.3 


1/512 


0.03 


0.02 


5.8 


4.6 



TABLE IV: Energy level splittings with and without exchange 
and correlation. 

closer to 1. In the regime above 1/16 ML, the modified 
potential results in a change of less than 6% in dp, and 
less than 12% in <za- 

The filling fractions, /3\, fa and fa\, are relatively in- 
sensitive to this perturbation. For 1/4 ML doping den- 
sity as described above, and indeed for densities above 
1/16, the values are almost identical with and without 
exchange and correlation, although there is a marked dif- 
ference at low densities. Table III details values of fa , fa 
and fa\. 

Table IV shows the effects of exchange and correla- 
tion, and density on the various energy level splittings. 
As might be expected, the valley-splitting also shows lit- 
tle effect of the inclusion or exclusion of exchange and 
correlation for those doping densities where Vxc is less 
significant. For densities higher than 1/16 ML, the dif- 
ference in the relative predicted splitting is less than 6%. 
Of course, if instead one is interested in the absolute dif- 
ference in the predictions, then exchange and correlation 



Dopant density 


Ef 


Ef 


(ML) 


(meV) (meV, no XC) 


1 


474 


668 


2/3 


251 


414 


1/2 


151 


296 


1/3 


63.1 


185 


1/4 


24.6 


132 


1/8 


-22.5 


57.7 


1/16 


-35.7 


24.0 


1/32 


-35.5 


8.97 


1/64 


-30.6 


2.81 


1/128 


-24.8 


0.32 


1/256 


-19.3 


-0.38 


1/512 


-14.8 


-0.57 



TABLE V: Fermi energy values with and without exchange 
and correlation. Note that the asymptotic value of the poten- 
tial, V (oo), has been defined as the energy zero. 



lead to an increase in the splitting of more than 0.1 meV 
for all doping densities above 1/32 ML. 



The Ti-A energy gap behaves similarly, with a rela- 
tive difference of less than 2% for systems denser than 
1/16 ML. As does the r 1 -r 2 splitting, the energy gap in- 
creases with inclusion of exchange and correlation. This 
is unsurprising, since the deeper potential corresponds to 
stronger confinement and hence wider spacings between 
energy levels. 



Table V displays the calculated Fermi energy with and 
without exchange and correlation. It is of note that we 
observe Fermi energies greater than zero, corresponding 
to an unphysical "overfilling" , over a much larger range 
of densities when ignoring exchange and correlation, with 
all densities above 1/128 ML (inclusive) having positive 
Fermi energies. This effect is more pronounced for higher 
densities. 



We note that, despite the Fermi energies changing 
by several meV for the 1/16-1/4 ML models, and be- 
ing greater than zero when exchange and correlation are 
excluded, all other results are largely unaffected by the 
inclusion or exclusion of exchange and correlation in the 
calculation. This is in line with the observations made 
in Ref. 26, where the Fermi level was by far the measure 
most sensitive to exclusion of exchange and correlation. 
We therefore observe that the overfilling appears to have 
little effect on these results - again in line with the dis- 
cussion regarding the populations of higher bands in Ref. 
26. We may then also consider it to have a similarly small 
effect on the higher densities, which also exhibit change 
due to the inclusion of exchange and correlation. 



While to first order, the inclusion of exchange and cor- 
relation appears to have little effect on our main results, 
we would like to emphasize its contribution to second- 
order effects such as the Fermi energy. We have noted 
that for several considered densities, including exchange 
and correlation reduces or even eliminates overfilling. It 
can easily be imagined that, when connected in a phys- 
ical device, significant overfilling would lead to breaking 
charge neutrality as the high-energy electrons are ener- 
getically free to vanish into the leads. As charge neu- 
trality is a central assumption in the derivation of our 
model, this is more important than perhaps it first ap- 
pears. We therefore recommend the inclusion of correla- 
tion and exchange in any EMT model of this type, espe- 
cially if larger-scale device modeling is to be undertaken. 
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